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1.  Introduction 

Outdoor  acoustic  propagation  is  strongly  affected  by 
the  covering  over  the  ground,  such  as  grass,  snow,  forest, 
or  asphalt.  To  extend  the  usefulness  of  acoustic  sensors 
to  urban  environments,  we  need  to  understand  and 
accurately  simulate  how  acoustic  waves  propagate  in  the 
presence  of  many  reflective  surfaces  offered  by  buildings 
and  man-made  ground  surfaces.  Work  by  others 
modeling  the  physical  mechanisms  of  how  acoustic 
energy  is  absorbed  by  various  materials  is  now  coded  into 
an  acoustic  propagation  code  and  used  to  simulate  how 
sound  propagates  over  various  types  of  terrain  and 
infrastructure  covered  by  a variety  of  materials. 

2.  Overview 

In  this  paper  we  review  the  theoretical  basis  for  the 
PSTOP3D  acoustic  propagation  code  and  the  recently- 
added  features  for  acoustic  source  modeling  and 
propagation  over  porous  surfaces.  We  then  present 
results  from  verification  efforts  for  free-field  propagation, 
propagation  over  porous  surfaces,  and  comparison  with 
experimental  data  for  reflections  from  a wall.  Finally,  we 
review  recent  efforts  to  extend  the  modeling  capabilities 
in  the  code,  as  well  as  continuing  modeling  efforts. 

3.  Description  of  PSTOP3D 

Our  computations  are  performed  using  a three- 
dimensional  (3D)  finite-  difference  time-domain  (FDTD) 
code.  The  code  suite  is  portable  and  is  operational  on 
several  Department  of  Defense  (DoD)  High  Performance 
Computing  Modernization  Program  (HPCMP)  machines, 
and  is  used  to  create  high  fidelity  data  sets  of  acoustic  or 
seismic  waves  propagating  over  realistic  terrain,  and 
includes  the  ability  to  model  buildings  and  infrastructure. 
For  acoustics,  the  equations  of  motion  and  pressure  are 
expressed  in  a rectangular  variable-grid  system  and 
discretized  with  second-order  differences.  The 


discretization  scheme  is  staggered  in  time  and  space, 
allowing  the  use  of  an  integration  scheme  known  as  the 
leapfrog  approach1 1].  Properties  at  boundaries  between 
layers  are  averaged  using  a scheme  developed  and 
demonstrated  by  Schroder  and  Scott121  to  enhance  stability 
across  the  boundaries  of  high-contrast  materials. 

The  size  (total  number  of  grid  nodes)  of  the  models 
that  can  be  run  is  proportional  to  the  number  of  processors 
and  memory  available  to  each  processor.  The  physical 
size  of  the  model  is  proportional  to  the  desired  bandwidth 
of  the  acoustic  response  to  be  simulated.  The  run  time  is 
proportional  to  the  length  of  time  of  the  simulation.  An 
example  which  was  run  on  Sapphire  (the  Engineer 
Research  and  Development  Major  Shared  Resource 
Center’s  [ERDC  MSRC’s]  XT3,  2007  upgrade,  with  two 
cores  per  processor  with  two  GB  ram  per  core)  is  shown 
in  Table  1 . Resources  required  by  PSTOP3D  to  run  much 
larger  acoustic  models  are  listed  in  Parker,  et  al.[31 
Example  resource  requirements  are  shown  in  Table  1. 

Table  1.  PSTOP3D  Resource  Requirements 


Parameter 

Job  1 

Job  2 

Total  nodes  (millions) 

1 

257 

Grid  spacing  (m) 

0.25 

0.25 

Physical  size  (m) 

25x25*25 

288x88x88 

Time  step  (s) 

2500 

1800 

Number  of  cores 

4 

CM 

CO 

Wall  time  (s) 

166 

216,000 

There  are  other  methods  for  calculating  the  acoustic 
response  over  large  scale  terrain.  These  include  ray- 
tracing methods141,  finite-element  modeling  approach151, 
and  the  Fast  Field  Program161.  The  advantage  of  the 
FDTD  approach  is  it  is  very  easy  to  mesh  the  problem 
domain,  and  algorithm  is  very  simple.  For  simulating 
acoustic  propagation  requires  only  the  implementation  of 
these  two  relations  (expanded  to  three-dimensions  from 
two-dimension  description  in  Wilson  and  Liul?1). 
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~ = V v + KeQ  (1) 

tit 

^-  = -—{Vp  + av)  (2) 

at  Pe 

where  p is  the  acoustic  pressure,  v is  the  vector  of  particle 
velocity  components,  t is  time,  Ke  is  the  effective  bulk 
modulus  in  the  material,  V is  the  divergence  operator,  pt 
is  the  effective  density  of  the  medium,  and  a is  the  static 
flow  resistivity  of  the  material  (which  is  by  definition 
zero  for  air,  and  nonzero  for  the  porous  materials  used  for 
ground  and  structures — this  is  described  further  in 
Section  5).  In  this  form  of  the  equations,  the  only  input 
term,  Q,  represents  a source  idealized  as  a radially 
oscillating  sphere®.  The  units  of  the  input  are  change  in 
volume  of  the  sphere  per  unit  time  per  unit  volume.  We 
will  say  more  about  the  source  term  in  the  next  section. 

The  disadvantage  of  the  FDTD  approach  is  the 
program  is  computationally  expensive,  in  that  it  requires 
node  spacing  of  about  10  nodes  per  wavelength,  and  the 
response  at  every  single  node  point  in  the  domain  is 
calculated  before  moving  to  the  next  time  step.  As  a side 
note,  the  finite-element  approach  being  developed  as 
described  in  Reference  5 uses  fewer  elements  per 
wavelength  and  may  provide  a more  efficient  approach. 

The  PSTOP3D  program  has  been  in  use  for  several 
years  at  ERDC,  and  has  been  used  to  study  acoustic  and 
seismic  propagation  over  realistic  terrain191,  tracking  of 
vehicles1®,  detection  of  tunnels  and  underground 
facilities,  and  simulation  of  the  effectiveness  of 
sensors11 'l  The  program  has  been  used  for  private-sector 
customers  under  a Cooperative  Research  and 
Development  Agreement  (CRADA)  and  distributed — 
with  accompanying  training — to  US  Army  Research  and 
Development  Engineering  Center  (ARDEC)  for  use  in  the 
Close  Combat  Support  Systems.  The  program  is 
undergoing  continuously  being  improved  in  terms  of  ease 
of  use,  verification,  validation,  documentation,  and 
incorporating  the  physics  necessary  to  describe  ever-more 
detailed  propagation  problems. 

4.  Acoustic  Source  Modeling 

The  only  input  term  to  Eqs.  (1  and  2)  is  the  source 
term  Q.  It  is  an  idealized  source  in  that  it  represents  a 
sphere  oscillating  about  some  nominal  radius  with  some 
velocity,  affecting  the  surrounding  medium  equally  in  all 
directions®.  The  radial  vibration  of  the  sphere  serves  to 
alternately  compress  and  rarefy  the  medium,  and  this 
change  in  density  corresponds  to  a change  in  pressure 
through  the  ideal  gas  law.  However,  the  volumetric 
change  per  unit  time  per  unit  volume  for  various  sources 
is  not  generally  known.  We  typically  measure  changes  in 


pressure  using  microphones,  but  the  measured  changes  in 
pressure  are  an  effect,  not  a cause,  of  a physical  process 
which  creates  sound.  Treating  measured  pressure 
fluctuation  as  an  input  to  the  equations  doesn’t  work  from 
the  point  of  view  of  the  units  of  the  equation,  or  from 
considerations  of  causality. 

Nonetheless,  we  understand  most  sources  in  terms  of 
the  pressure  fluctuations  they  cause.  Pressure  measured 
at  a distance  from  a source  can  be  modeled  as  a 
volumetric  source  through  the  relation®, 

Qs{t)~—\p{t)dt  (3) 

P J 


where  (9s(t)  is  the  required  change  in  volume  of  the  source 
per  unit  time  which  would  cause  the  pressure  fluctuations 
p(t),  and  r is  the  distance  of  the  pressure  measurement 
from  the  source.  Equation  (3)  has  the  units  m3/s,  and  is 
the  total  volumetric  source.  This  equation  provides  a 
means  of  converting  pressure  measured  at  a distance  to  an 
ideal  source,  limited  by  the  following  two  assumptions: 
1)  the  source  is  modeled  as  acting  at  a point  a distance  r 
away  from  the  measurement,  and  2)  the  pressure 
measured  at  the  distance  r is  assumed  to  be  constant  over 
a spherical  surface  of  radius  r. 

Currently,  PSTOP3D  can  accommodate  multiple 
sources,  but  treats  acoustic  sources  as  acting  at  a single 
grid  point.  For  the  large-domain  problems  the  grid 
spacing  can  be  from  0.1  to  0.3  m,  and  the  point  source 
assumption  is  sufficiently  accurate  for  many  sources  of 
interest  (e.g.,  small-arms  fire,  small  explosions).  So  the 
first  assumption  is  not  a severe  limitation  given  our 
spatial  scales  and  limited  frequency  range  for  large  scale 
simulations.  Furthermore,  distributing  the  acoustic 
sources  over  a domain  which  reflects  the  size  of  the 
source  would  be  straightforward. 

The  second  assumption  is  more  of  a limitation. 
Many  sources  are  measured  in  the  free  field,  with  the 
microphone  often  1 m from  the  source.  But  most  sources 
have  some  directivity,  and  the  source  power  is 
characterized  in  general  terms.  The  FDTD  approach  can 
accommodate  sources  with  directivity,  but  the  necessary 
free-field  input  data  is  not  often  available  for  large 
physical  sources  such  as  vehicles  operating  over  terrain. 

The  ideal  volumetric  source  Q{i)  which  causes  the 
measured  pressure  at  the  distance  r must  have  the  units 
1/s  to  satisfy  Eq.  (1),  and  this  means  dividing  the  total 
required  change  in  volume  of  the  source  per  unit  time  by 
the  volume  over  which  the  source  is  defined  to  get  the 
change  of  volume  per  unit  time  per  unit  volume.  In  the 
FDTD  approach,  for  a source  acting  at  a single  grid  point, 
this  means  dividing  by  the  unit  volume  at  that  grid  point: 


0(0= 


0,(0 

dx-dy-dz 


(4) 
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where  dx,  dy,  and  dz  are  the  dimensions  of  the  cubo-id 
quadrilateral  hexahedron  (if  dx,  dy,  and  dz  are  equal,  this 
is  simply  a cube). 

5.  Porous  Surface  Modeling 


The  term  av  in  Eq.  (2)  represents  one  effect  of 
porosity.  The  effect  of  porous  surfaces  on  acoustic 
propagation  has  been  performed  for  2D  simulations, 
(Wilson,  et  al.[12]),  based  on  the  work  by 
Attenborough1'3'141.  We  simply  extend  this  to  3D  and 
investigate  constraints  imposed  by  considerations  of 
stability. 

The  physics  of  acoustic  propagation  in  the  free  field 
is  reflected  in  Eqs.  (1  and  2),  if  term  cis  zero.  To  include 
the  effects  of  porous  materials  on  acoustic  propagation, 
the  following  parameters  are  used:  the  void  fraction  Q 
(also  called  the  material  porosity,  and  is  the  volume  of  the 
pores  per  unit  volume,  m3/m3),  static  flow  resistivity  <r 
(which  reflects  the  relative  size  of  the  pores,  and  which 
has  the  units  Pa/m2/s),  and  tortuosity,  q,  (no  units)  which 
reflects  the  geometry  of  the  pores.  These  parameters  were 
chosen  as  the  simplest  and  most  effective  implementation 
of  the  extensive  research  on  the  effect  of  porosity  on 
acoustic  propagation.  Wilson  and  Liu[71  provide  a 
convenient  tabulation  of  these  parameters  tortuosity, 
porosity,  and  flow  resistivity  for  asphalt,  snow,  sand, 
grass,  and  forest. 

We  have  implemented  equations  for  acoustic 
propagation  for  porous  materials  based  on  the  following 
assumptions:  (1)  the  pores  are  cylindrical,  with  the  axial 
axis  of  the  cylinder  normal  to  the  surface  of  the  material, 
(2)  our  analysis  is  limited  to  relatively  low  frequencies 
(where  the  acoustic  wavelengths  are  much  larger  than  the 
pore  dimensions),  and  (3)  there  is  negligible  temperature 
change  due  to  the  acoustic  propagation,  so  isothermal 
relations  apply). 

To  go  from  tortuosity  and  porosity  to  the  effective 
bulk  modulus  and  effective  density,  the  following 
relations  are  used  (limited  to  the  assumptions  above  and 
modified  slightly  from  Reference  7): 


Ke 


(5) 


4 pq1 

3 Q 


(6) 


where  yis  the  ratio  of  specific  heats  of  fluid  portion  of  the 
porous  material,  c is  the  speed  of  sound  in  the  fluid 
portion  of  the  porous  material,  and  p is  the  density  of  the 
fluid  portion  of  the  porous  material. 

The  effect  of  porous  materials  are  implemented  in  the 
preprocessor  code  portion  of  PSTOP3D  in  such  a way 
that,  during  the  high-performance  computing  processing, 


Eqs.  (1  and  2)  are  implemented  across  the  whole 
computational  domain.  This  is  possible  because  the 
equations  for  porous  materials  are  the  same  as  single-state 
(e.g.,  all  gas,  all  liquid,  or  all  solid)  if,  in  the  preprocessor, 

the  effective  terms  K and  pe  are  set  equivalent  to  k and 
p,  respectively,  and  the  static  flow  resistivity  is  set  to 
zero. 

Our  integration  scheme  is  based  on  the  necessary  and 
sufficient  condition  for  the  size  of  the  time  step  during 
integration  to  be  bounded  by  the  Courant  condition. 
However,  when  highly  reflective  porous  materials  were 
implemented  (e.g.,  asphalt),  the  simulations  went 
unstable,  and  the  instability  started  where  the  propagating 
wave  first  encountered  the  porous  material.  (Note,  as 
mentioned  previously,  we  have  implemented  the  results 
from  Reference  2 to  maximize  stability  at  an  elastic 
interface).  The  results  are  explained  if  we  examine  the 
difference  equation  form  of  Eq.  (2)  used  to  update  the 
particle  velocity  from  the  old  time  step  to  the  next  time 
step: 


' A/  ' 
1 cr 

Pc  J 


A? 

'aid VP 

Pe 


(7) 


where  At  is  the  time  step(s).  If  we  consider  Eq.  (7)  as  a 
finite-difference  form  of  a first  order  equation  in  v (by 
taking  Vp  to  be  constant  over  the  time  step),  we  recognize 


that  for  stability,  the  term 


Ar 

l a 

V P. 


must  be  greater  than 


zero. 

We  can  meet  this  condition  by  limiting  At,  but 
materials  such  as  asphalt  have  a static  flow  resistivity 
equal  to  3e7  Pa/m2/s.  Setting  At  to  meet  this  condition 
would  increase  the  computational  cost  by  orders  of 
magnitude. 

We  also  investigated  using  central  differences  to 
approximate  Eqs.  (1  and  2);  however,  this  approach 
requires  more  memory  (to  store  the  intermediate  values  of 
the  fields  for  pressure  and  velocity;  and  this  storage 
requirement  triples  in  the  seismic  propagation 
simulations),  and,  although  more  stable  than  the  Euler 
method  implemented  in  PSTOP3D  (in  that  it  typically  can 
use  a larger  time  step  and  remain  stable171),  there  was  still 
an  onerous  computational  burden  to  implement  the  time 
step  required  for  a stable  solution.  Please  note  that  a 
specific  condition  for  stability  using  the  central  difference 
approach  was  not  derived. 

The  solution  implemented  that  met  both  the 
resources-  based  requirement  of  using  the  Courant 
condition  to  set  the  time  step,  and  the  accuracy 
requirement  for  simulating  over  long  distances,  was  to 
artificially  cap  the  value  of  the  static  flow  resistivity  to 
satisfy, 
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(8) 


7.  Propagation  over  Porous  Layers 


Pe 

The  results  show  good  agreement  for  asphalt,  which 
is  the  most  reflective  of  materials  we  implemented,  and 
which  is  most  significantly  impacted  by  implementing 
Eq.  (8).  Equation  (8)  is  the  relationship  that  was  used  for 
simulations  shown  in  References  3 and  1 1 . 

6.  Ideal  Free-field  Propagation 

The  purpose  of  this  verification  is  to  ensure  that 
simulations  performed  using  PSTOP3D  for  propagation  in 
the  free-field  agree  with  theory  over  distances  up  to  200 
m,  and  that  our  modeling  approach  for  modeling  sources 
based  on  a measured  pressure  yields  the  exact  measured 
pressure  at  the  same  distance  from  the  source. 

A small  simulation  for  examining  the  source  used  a 
40  Hz  harmonic  pulse  source  located  in  the  geometric 
center  of  a model  100  nodes  on  each  side,  with  an 
absorbing  layer  of  10  nodes  on  all  six  sides,  using  a 
cosine-based  grid  stretching  formula  with  an  amplitude  of 
4,  and  the  model  being  completely  composed  of  air  with 
an  acoustic  wave  speed  of  347  m/s  and  a density  of  1 .2 
kg/m3.  To  show  that  the  results  are  independent  of  grid 
spacing,  two  different  grid  spacings  are  used,  0.25  and  0. 1 
m,  with  corresponding  integration  time  steps  of  0.414  and 
0.165  ms,  respectively.  The  results  are  shown  in  Figure  1, 
and  the  acoustic  source  term  generates  the  required 
pressure  amplitude  of  1 Pa  at  1 m from  the  source. 

A larger  simulation  was  run  using  a grid  of  1280  x 
448  x 448  nodes,  with  a grid  spacing  or  0.25  m,  with  a 
harmonic  source  located  44  m away  from  five  boundaries. 
The  air  properties  and  the  absorbing  layer  are  the  same  as 
the  previous  example.  The  sinusoidal  harmonic  is 
analytically  integrated  to  give  a cosine  function,  offset  to 
start  at  0 m3/m3/s,  and  filtered  to  avoid  starting  transients 
for  the  higher-frequency  inputs.  The  filtered  input  signal 
is  shown  in  Figure  2 for  the  40  Hz  harmonic  input,  and 
the  results  compared  to  theory  for  transmission  loss  in  air 
are  shown  in  Figures  3 and  4,  for  harmonic  sources  at  40 
and  100  Hz,  respectively.  Note  that  the  definition  for 
sound  transmission  loss  is, . 

TL  = 201og  1 0[p(x)/p(x=  1 m)]  (9) 

where  p(x)  is  the  root-mean-square  of  the  pressure 
amplitude.  The  theory  for  pressure  amplitude  as  a 
function  of  distance  from  the  source  is, 

co) 

The  results  show  excellent  agreement,  with  errors 
building  up  to  less  than  1 dB  200  m from  the  source  for 
the  100  Hz  harmonic  input. 


Initial  results  comparing  propagation  over  porous 
layers  is  shown  in  Figure  5.  The  parameters  of  the 
problem  are  shown  in  Table  2.  The  Fast  Field  Program 
has  been  extensively  verified  and  tested  for  propagation 
over  porous  materials.  Note  that  the  model  capped  the 
static  flow  resistivity  in  a slightly  different  way. 

Table  2.  Problem  parameters  for  verifying  accuracy  of 
propagation  due  to  asphalt 


Parameter 

Value 

Grid 

1,400  x 1,400x580 

Spacing  (m) 

0.25 

Domain  (unstretched,  m) 

350  x 350  x 145 

Absorbing  region 

20  elements,  scale 
factor  = 4 

Time  step  (ms) 

0.414 

Pressure  @ 1 m from  source 
(Pa) 

1.0 

Source  location  - three 
different  heights  (m) 

x=  150.2,  y=  177.2, 
z=[24.2;  27.9;  42.9] 

Frequency  (Hz) 

50, 100  Hz 

Receiver  locations  - three 
different  heights,  extended 
series  in  x (m) 

150.2  + 0.5  m*400 
points,  177.2,  [24.2; 
30.2;  45.2] 

Ground  height  (m) 

21.9 

8.  Response  Due  to  an  Explosion 

Based  on  an  experiment  by  Albert  & Liu[l5],  we  have 
an  analytical  model  which  closely  approximates  the 
pressure  from  a high-order  explosion  of  a single  block  of 
C4.  The  equation  for  the  pulse  is, 

p{t)  = |l-16[l50-/-0.25]2  je~10(150',)/3  (11) 

The  results  are  shown  in  Figure  6. 

9.  Summary 

By  modeling  pressure  using  an  idealized  oscillating 
sphere,  we  can  reproduce  the  exact  pressure  measured  in 
the  free  field.  Due  to  the  computational  costs  required  to 
meet  the  stability  constraints  we  are  forced  to  artificially 
limit  the  static  flow  resistivity,  however,  comparison  to 
the  Fast  Field  Program  for  asphalt  (which  has  the  highest 
static  flow  resistivity  in  the  group  of  materials  for  which 
we  have  values),  we  see  the  predicted  reflectivity  has  very 
little  error  when  compared  to  an  ideal  rigid  (i.e., 
completely  reflective)  surface. 
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Further  work  will  extend  the  verification  effort  to 
other  materials,  validate  the  code  by  comparing  to 
experimental  data,  and  develop  a method  to  model  an 
ideal  source  based  on  pressure  measurements  which  are 
not  completely  in  the  free-field. 
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Figure  1.  Response  to  40  Hz  harmonic  pulse 
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Figure  2.  Input  Q(t)  = 100  Hz  sine  wave,  after 
integration  and  low-pass  filtering 
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Figure  3.  Comparison  of  PSTOP3D  to  Fast  Field 
Program,  40  Hz  harmonic  input 
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Figure  4.  Comparison  of  PSTOP3D  to  Fast  Field 
Program,  100  Hz  harmonic  input 


Figure  5.  Reflectivity  of  asphalt  calculated  by 
PSTOP3D  and  Fast  Field  Program,  compared  to 
propagation  over  an  ideal  rigid  surface 
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Figure  6.  Comparison  of  output  from  PSTOP3D  with 
original  pulse.  (Original  pulse  shifted  to  account  for 
delay  introduced  by  filtering). 
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